Geospatial Analysis using OneMap Data

Published

February 24, 2024

Data Analysis

pacman::p_load(arrow, lubridate, maptools, sp, sf, raster, spatstat, tmap, classInt, viridis, tidyverse, spNetwork)

Reading shapefiles that were exported from our qgis geopackage:

pg_wdl <- st_read(dsn = "data/geospatial", 
                 layer = "pg_wdl")
Reading layer `pg_wdl' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 3938 features and 12 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34031.82 ymin: 40392.88 xmax: 38881.06 ymax: 44812.07
Projected CRS: SVY21 / Singapore TM
amk_wdl <- st_read(dsn = "data/geospatial", 
                 layer = "wdl_amk")
Reading layer `wdl_amk' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 6855 features and 12 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 26152.2 ymin: 37505 xmax: 31404.37 ymax: 42104.04
Projected CRS: SVY21 / Singapore TM
dl_pg <- st_read(dsn = "data/geospatial", 
                 layer = "pg_dl")
Reading layer `pg_dl' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 40 features and 12 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 35253.61 ymin: 42595.38 xmax: 37130.03 ymax: 43050.74
Projected CRS: SVY21 / Singapore TM
dl_amk <- st_read(dsn = "data/geospatial", 
                 layer = "dl_amk")
Reading layer `dl_amk' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 13 features and 12 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 28451.72 ymin: 38451.37 xmax: 30392.43 ymax: 39764.13
Projected CRS: SVY21 / Singapore TM

library(ggplot2)
#install.packages("ggmap")
library(ggmap)
#dl_map2 <- osm_basemap +
#  tm_shape(dl_class) +
#  tm_lines(line_colors)
#  tmap_options(check.and.fix = TRUE)

# View the map
#dl_map2
target_buildings <- st_read(dsn = "data/geospatial", 
                 layer = "target_buildings")
Reading layer `target_buildings' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 9554 features and 7 fields (with 1 geometry empty)
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 26581.25 ymin: 33134.58 xmax: 44722.63 ymax: 44778.11
Projected CRS: SVY21 / Singapore TM
Coordinate Reference System:
  User input: SVY21 / Singapore TM 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]
tb_cat <- target_buildings["Categorize"]

Reading layer `MPSZ-2019' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial/MPSZ' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
mpsz3414 <- st_transform(mpsz, 3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

pa_mpsz <- mpsz3414["PLN_AREA_N"]
mpsz2 <- as_Spatial(pa_mpsz)
pg = mpsz2[mpsz2@data$PLN_AREA_N == "PUNGGOL",]
amk = mpsz2[mpsz2@data$PLN_AREA_N == "ANG MO KIO",]

pg_sp = as(pg, "SpatialPolygons")
amk_sp = as(amk, "SpatialPolygons")
pg_owin = as(pg_sp, "owin")
amk_owin = as(amk_sp, "owin")
edl_pg_points_sf <- st_as_sf(edl_pg_points)
edl_coords <- st_coordinates(edl_pg_points_sf)

# Create an owin object from the coordinates
edl_pg_owin <- owin(xrange = range(edl_coords[,1]), yrange = range(edl_coords[,2]))
edl_amk_points_sf <- st_as_sf(edl_amk_points)
amk_edl_coords <- st_coordinates(edl_amk_points_sf)

# Create an owin object from the coordinates
edl_amk_owin <- owin(xrange = range(amk_edl_coords[,1]), yrange = range(amk_edl_coords[,2]))

edl_pg_ppp = edl_pg_ppp[pg_owin]
edl_a_ppp = edl_amk_ppp[amk_owin]

class(pg)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
pg_sf <- st_as_sf(pg)
amk_sf <- st_as_sf(amk)
pg_buildings <- st_intersection(pg_sf, tb_cat)
amk_buildings <- st_intersection(amk_sf, tb_cat)

library(ggplot2)
#install.packages("ggmap")
library(ggmap)
landuse_pg <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_target_landuse_PG")
Reading layer `FYP_target_landuse_PG' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 567 features and 22 fields
Geometry type: MULTIPOLYGON
Dimension:     XY, XYZ
Bounding box:  xmin: 31956.1 ymin: 40850.4 xmax: 38889.96 ymax: 44808.79
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
landuse_amk <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_target_landuse_AMK")
Reading layer `FYP_target_landuse_AMK' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 5018 features and 22 fields
Geometry type: MULTIPOLYGON
Dimension:     XY, XYZ
Bounding box:  xmin: 26147.34 ymin: 37481.8 xmax: 31185.22 ymax: 42180.82
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
landuse_pg_desc <- (landuse_pg["descriptio"])
plot(landuse_amk["descriptio"])

plot(landuse_pg_desc)

landuse_pg_sf <- st_as_sf(landuse_pg)
landuse_amk_sf <- st_as_sf(landuse_amk)
landuse_pg_sf <- st_make_valid(landuse_pg_sf)
valid_pg <- st_make_valid(pg_sf)
valid_amk <- st_make_valid(amk_sf)
pg_landuse <- st_intersection(valid_pg, landuse_pg_sf)
amk_landuse <- st_intersection(valid_amk, landuse_amk_sf)
parks_pg <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_walkingtime_parks_PG")
Reading layer `FYP_walkingtime_parks_PG' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 20 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 34876.04 ymin: 41522.67 xmax: 37386.37 ymax: 43523.24
Projected CRS: SVY21 / Singapore TM
parks_amk <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_walkingtime_parks_AMK")
Reading layer `FYP_walkingtime_parks_AMK' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 11 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 27351.72 ymin: 39057.31 xmax: 30492.3 ymax: 41605.77
Projected CRS: SVY21 / Singapore TM
healthcare_pg <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_walkingtime_healthcare_PG")
Reading layer `FYP_walkingtime_healthcare_PG' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 21 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 34987.76 ymin: 41965.05 xmax: 38171.03 ymax: 44135
Projected CRS: SVY21 / Singapore TM
healthcare_amk <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_walkingtime_healthcare_AMK_1")
Reading layer `FYP_walkingtime_healthcare_AMK_1' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 16 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 27176.45 ymin: 38409.92 xmax: 31090.67 ymax: 39967.57
Projected CRS: SVY21 / Singapore TM
food_pg <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_walkingtime_foodplaces_PG")
Reading layer `FYP_walkingtime_foodplaces_PG' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 12 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 34861.09 ymin: 41878.28 xmax: 37764.21 ymax: 43760.28
Projected CRS: SVY21 / Singapore TM
food_amk <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_walkingtime_foodplaces_AMK")
Reading layer `FYP_walkingtime_foodplaces_AMK' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 14 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 28389.21 ymin: 38453.36 xmax: 30971.45 ymax: 40817.5
Projected CRS: SVY21 / Singapore TM
smarket_pg <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_walkingtime_smarket_PG")
Reading layer `FYP_walkingtime_smarket_PG' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 16 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 34958.94 ymin: 41808.02 xmax: 37583.44 ymax: 44161.24
Projected CRS: SVY21 / Singapore TM
smarket_amk <- st_read(dsn = "data/geospatial", 
                 layer = "FYP_walkingtime_smarket_AMK")
Reading layer `FYP_walkingtime_smarket_AMK' from data source 
  `/Users/shanellelam/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/School/Year 3/chanelelele/fyp_newnetli_shan/Analysis/Geographic Analysis/data/geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 12 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 27166.86 ymin: 38567.79 xmax: 30312.63 ymax: 40817.5
Projected CRS: SVY21 / Singapore TM
plot(food_pg)

parks_pg_sf <- st_as_sf(parks_pg)
parks_amk_sf <- st_as_sf(parks_amk)
healthcare_pg_sf <- st_as_sf(healthcare_pg)
healthcare_amk_sf <- st_as_sf(healthcare_amk)
food_pg_sf <- st_as_sf(food_pg)
food_amk_sf <- st_as_sf(food_amk)
smarket_pg_sf <- st_as_sf(smarket_pg)
smarket_amk_sf <- st_as_sf(smarket_amk)
parks_pg_sf <- st_make_valid(parks_pg_sf)
parks_amk_sf <- st_make_valid(parks_amk_sf)
healthcare_pg_sf <- st_make_valid(healthcare_pg_sf)
healthcare_amk_sf <- st_make_valid(healthcare_amk_sf)
food_pg_sf <- st_make_valid(food_pg_sf)
food_amk_sf <- st_make_valid(food_amk_sf)
smarket_pg_sf <- st_make_valid(smarket_pg_sf)
smarket_amk_sf <- st_make_valid(smarket_amk_sf)
pg_parks <- st_intersection(valid_pg, parks_pg_sf)
pg_healthcare <- st_intersection(valid_pg, healthcare_pg_sf)
pg_food <- st_intersection(valid_pg, food_pg_sf)
pg_smarket <- st_intersection(valid_pg, smarket_pg_sf)

amk_parks <- st_intersection(valid_amk, parks_amk_sf)
amk_healthcare <- st_intersection(valid_amk, healthcare_amk_sf)
amk_food <- st_intersection(valid_amk, food_amk_sf)
amk_smarket <- st_intersection(valid_amk, smarket_amk_sf)

Punggol Buildings Chart

pg_b_cat <- pg_buildings["Categorize"]
plot(pg_b_cat, key.width = lcm(6.32))

Punggol Buildings Chart on OSM Layer

tmap_mode('view')
# Define basemaps
osm_basemap <- tm_basemap(server = "OpenStreetMap.HOT")
imagery_basemap <- tm_basemap(server = "Esri.WorldImagery")

# Plot pg_b_cat over OpenStreetMap basemap with color categorization and legend
pg_b_cat_map2 <- tm_shape(pg_b_cat) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "Categorize", title = "Building Category") +  # Color categorization
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
pg_b_cat_map2

Ang Mo Kio Buildings Chart

amk_b_cat <- amk_buildings["Categorize"]
plot(amk_b_cat, key.width = lcm(6.32))

Ang Mo Kio Buildings Chart on OSM Layer

# Plot pg_b_cat over OpenStreetMap basemap with color categorization and legend
amk_b_cat_map <- tm_shape(amk_b_cat) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "Categorize", title = "Building Category") +  # Color categorization
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
amk_b_cat_map

Punggol Walkway Lines + Desired Lines on OSM Layer

tmap_mode('view')

# Define basemaps
osm_basemap <- tm_basemap(server = "OpenStreetMap.HOT")
imagery_basemap <- tm_basemap(server = "Esri.WorldImagery")

# Plot dl_class over OpenStreetMap basemap
dl_map <- osm_basemap +
  tm_shape(dl_class) +
  tm_lines(col = "red", lwd = 2) +  # Adjust line color and width as needed
  tmap_options(check.and.fix = TRUE)

# View the map
dl_map

Ang Mo Kio Walkway Lines + Desired Lines on OSM Layer

a_dl_class <- amk_wdl["desired_li"]
tmap_mode('view')

# Define basemaps
osm_basemap <- tm_basemap(server = "OpenStreetMap.HOT")
imagery_basemap <- tm_basemap(server = "Esri.WorldImagery")

# Plot dl_class over OpenStreetMap basemap
dl_map <- osm_basemap +
  tm_shape(a_dl_class) +
  tm_lines(col = "red", lwd = 2) +  # Adjust line color and width as needed
  tmap_options(check.and.fix = TRUE)

# View the map
dl_map

Punggol Desired Lines vs Walkway Lines

# Define colors based on the value of the "fclass" column
line_colors <- ifelse(pg_wdl$fclass == "dl", "red", "green")

dl_class <- pg_wdl["desired_li"]

# Plot with custom colors
plot(dl_class, col = line_colors, main = "Desired Lines in Punggol")

Ang Mo Kio Desired Lines vs Walkway Lines

# Define colors based on the value of the "fclass" column
line_colors <- ifelse(amk_wdl$fclass == "dl", "red", "green")

a_dl_class <- amk_wdl["desired_li"]

# Plot with custom colors
plot(a_dl_class, col = line_colors, main = "Desired Lines in Ang Mo Kio")

Punggol Desired Lines on OSM Layer

tmap_mode('view')

# Define basemaps
osm_basemap <- tm_basemap(server = "OpenStreetMap.HOT")
imagery_basemap <- tm_basemap(server = "Esri.WorldImagery")

# Plot dl_class over OpenStreetMap basemap
dl_map <- osm_basemap +
  tm_shape(dl_pg) +
  tm_lines(col = "red", lwd = 2) +  # Adjust line color and width as needed
  tmap_options(check.and.fix = TRUE)

# View the map
dl_map

Ang Mo Kio Desired Lines on OSM Layer

tmap_mode('view')

# Define basemaps
osm_basemap <- tm_basemap(server = "OpenStreetMap.HOT")
imagery_basemap <- tm_basemap(server = "Esri.WorldImagery")

# Plot dl_class over OpenStreetMap basemap
dl_map_amk <- osm_basemap +
  tm_shape(dl_amk) +
  tm_lines(col = "red", lwd = 2) +  # Adjust line color and width as needed
  tmap_options(check.and.fix = TRUE)

# View the map
dl_map_amk

Punggol Kernel Density Estimation Graphs

# Set up a 2x2 layout with reduced margins
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))

# Chart 1
plot(density(edl_pg_ppp,
             sigma = bw.CvL,
             edge = TRUE,
             kernel = "gaussian"),
     main = "Punggol Desired Lines CvL")

# Chart 2
plot(density(edl_pg_ppp,
             sigma = bw.scott,
             edge = TRUE,
             kernel = "gaussian"),
     main = "Punggol Desired Lines scott")

# Chart 3
plot(density(edl_pg_ppp,
             sigma = bw.diggle,
             edge = TRUE,
             kernel = "gaussian"),
     main = "Punggol Desired Lines diggle")

# Chart 4
plot(density(edl_pg_ppp,
             sigma = bw.ppl,
             edge = TRUE,
             kernel = "gaussian"),
     main = "Punggol Desired Lines ppl")

Ang Mo Kio Kernel Density Estimation Graphs

# Set up a 2x2 layout with reduced margins
par(mfrow = c(2, 2), mar = c(2, 2, 2, 2))

# Chart 1
plot(density(edl_a_ppp,
             sigma = bw.CvL,
             edge = TRUE,
             kernel = "gaussian"),
     main = "Ang Mo Kio Desired Lines CvL")

# Chart 2
plot(density(edl_a_ppp,
             sigma = bw.scott,
             edge = TRUE,
             kernel = "gaussian"),
     main = "Ang Mo Kio Desired Lines scott")

# Chart 3
plot(density(edl_a_ppp,
             sigma = bw.diggle,
             edge = TRUE,
             kernel = "gaussian"),
     main = "Ang Mo Kio Desired Lines diggle")

# Chart 4
plot(density(edl_amk_ppp,
             sigma = bw.ppl,
             edge = TRUE,
             kernel = "gaussian"),
     main = "Ang Mo Kio Desired Lines ppl")

Punggol Landuse on OSM Layer

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
pg_landuse_map <- tm_shape(pg_landuse) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "Categorize", title = "Landuse Category", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
pg_landuse_map

Ang Mo Kio Landuse on OSM Layer

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
amk_landuse_map <- tm_shape(amk_landuse) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "Categorize", title = "Landuse Category", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
amk_landuse_map

Punggol Isotope Areas of Parks on OSM Layer (Key Amenities)

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
pg_parks_map <- tm_shape(pg_parks) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "cost_level", title = "Isotope Areas of Parks", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
pg_parks_map

Punggol Isotope Areas of Healthcare on OSM Layer (Key Amenities)

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
pg_healthcare_map <- tm_shape(pg_healthcare) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "cost_level", title = "Isotope Areas of Healthcare", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
pg_healthcare_map

Punggol Isotope Areas of Food Places on OSM Layer (Key Amenities)

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
pg_food_map <- tm_shape(pg_food) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "cost_level", title = "Isotope Areas of Food Places", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
pg_food_map

Punggol Isotope Areas of Supermarkets on OSM Layer (Key Amenities)

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
pg_smarket_map <- tm_shape(pg_smarket) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "cost_level", title = "Isotope Areas of Supermarkets", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
pg_smarket_map

Ang Mo Kio Isotope Areas of Parks on OSM Layer (Key Amenities)

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
amk_parks_map <- tm_shape(amk_parks) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "cost_level", title = "Isotope Areas of Parks", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
amk_parks_map

Ang Mo Kio Isotope Areas of Healthcare on OSM Layer (Key Amenities)

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
amk_healthcare_map <- tm_shape(amk_healthcare) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "cost_level", title = "Isotope Areas of Healthcare", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
amk_healthcare_map

Ang Mo Kio Isotope Areas of Supermarkets on OSM Layer (Key Amenities)

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
amk_smarket_map <- tm_shape(amk_smarket) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "cost_level", title = "Isotope Areas of Supermarkets", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
amk_smarket_map

Ang Mo Kio Isotope Areas of Food Places on OSM Layer (Key Amenities)

# Plot pg_landuse over OpenStreetMap basemap with color categorization, legend, and lower opacity
amk_food_map <- tm_shape(amk_food) +
  tm_borders() +  # Plot borders of the polygons
  tm_fill(col = "cost_level", title = "Isotope Areas of Food Places", alpha = 0.5) +  # Color categorization with lower opacity
  tm_layout(legend.show = TRUE) +  # Show legend
  osm_basemap  # Add basemap

# View the map
amk_food_map